<!--begin.rcode results='hide', echo=FALSE, message=FALSE
library(caret)
library(gbm)
library(pROC)
library(klaR)
library(kernlab)
theme1 <- caretTheme()
theme1$superpose.symbol$col = c(rgb(1, 0, 0, .4), rgb(0, 0, 1, .4), 
  rgb(0.3984375, 0.7578125,0.6445312, .6))
theme1$superpose.symbol$pch = c(15, 16, 17)
theme1$superpose.cex = 1
theme1$superpose.line$col = c(rgb(1, 0, 0, .9), rgb(0, 0, 1, .9), rgb(0.3984375, 0.7578125,0.6445312, .6))
theme1$superpose.line$lwd <- 2
theme1$superpose.line$lty = 1:3
theme1$plot.symbol$col = c(rgb(.2, .2, .2, .4))
theme1$plot.symbol$pch = 16
theme1$plot.cex = 1
theme1$plot.line$col = c(rgb(1, 0, 0, .7))
theme1$plot.line$lwd <- 2
theme1$plot.line$lty = 1

    hook_inline = knit_hooks$get('inline')
    knit_hooks$set(inline = function(x) {
      if (is.character(x)) highr::hi_html(x) else hook_inline(x)
      })
    opts_chunk$set(comment=NA, digits = 3)

library(ellipse)
library(rpart)
session <- paste(format(Sys.time(), "%a %b %d %Y"),
                 "using caret version",
                 packageDescription("caret")$Version,
                 "and",
                 R.Version()$version.string)

options(width = 150)

library(parallel)
library(doMC)
registerDoMC(cores=detectCores()-1)

    end.rcode--> 

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
  <!--
  Design by Free CSS Templates
http://www.freecsstemplates.org
Released for free under a Creative Commons Attribution 2.5 License

Name       : Emerald 
Description: A two-column, fixed-width design with dark color scheme.
Version    : 1.0
Released   : 20120902

-->
  <html xmlns="http://www.w3.org/1999/xhtml">
  <head>
  <meta name="keywords" content="" />
  <meta name="description" content="" />
  <meta http-equiv="content-type" content="text/html; charset=utf-8" />
  <title>Model Training and Tuning</title>
  <link href='http://fonts.googleapis.com/css?family=Abel' rel='stylesheet' type='text/css'>
  <link href="style.css" rel="stylesheet" type="text/css" media="screen" />
  </head>
  <body>
  <div id="wrapper">
  <div id="header-wrapper" class="container">
  <div id="header" class="container">
  <div id="logo">
  <h1><a href="#">Model Training and Tuning</a></h1>
</div>
  <!--
  <div id="menu">
  <ul>
  <li class="current_page_item"><a href="#">Homepage</a></li>
<li><a href="#">Blog</a></li>
<li><a href="#">Photos</a></li>
<li><a href="#">About</a></li>
<li><a href="#">Contact</a></li>
</ul>
  </div>
  -->
  </div>
  <div><img src="images/img03.png" width="1000" height="40" alt="" /></div>
  </div>
  <!-- end #header -->
<div id="page">
  <div id="content">

<h1>Contents</h1>  
<ul>
  <li><a href="#builtin">Model Training and Parameter Tuning</a></li>
  <ul>
     <li><a href="#example">An Example</a></li>
  </ul>   
  <li><a href="#tune">Basic Parameter Tuning</a></li>
  <li><a href="#repro">Notes on Reproducibility</a></li>  
  <li><a href="#custom">Customizing the Tuning Process</a></li> 
  <ul>
     <li><a href="#preproc">Pre-Processing Options</a></li>
     <li><a href="#grids">Alternate Tuning Grids</a></li>
     <li><a href="#plots">Plotting the Resampling Profile</a></li>
     <li><a href="#control">The <span class="mx funCall">trainControl</span> Function</a></li>
     <li><a href="#metrics">Alternate Performance Metrics</a></li>
     <li><a href="#final">Choosing the Final Model</a></li>
  </ul>   
  <li><a href="#pred">Extracting Predictions and Class Probabilities</a></li> 
  <li><a href="#resamp">Exploring and Comparing Resampling Distributions</a></li> 
  <ul>
     <li><a href="#within">Within-Model</a></li>
     <li><a href="#between">Between-Models</a></li>
  </ul>  
  <li><a href="#notune">Fitting Models Without Parameter Tuning</a></li> 
 </ul>   
      
<div id="builtin"></div> 
<h1>Model Training and Parameter Tuning</h1>

<p>
The <a href="http://cran.r-project.org/web/packages/caret/index.html"><strong>caret</strong></a> package has several functions that attempt to streamline the model building and evaluation process. <p>
</p> 
<p>
<p>
The <span class="mx funCall">train</span> function can be used to
  <ul>
  <li> evaluate, using resampling, the effect of model tuning parameters on performance
  <li> choose the "optimal" model across these parameters 
  <li> estimate model performance from a training set
  </ul>
</p> 

<p>
First, a specific model must be chosen. Currently,
<!--rinline I(length(unique(modelLookup()$model))) --> are available using
<a href="http://cran.r-project.org/web/packages/caret/index.html"><strong>caret</strong></a>; see <a href="modelList.html"><tt>train</tt> Model List</a> 
or <a href="bytag.html"><tt>train</tt> Models By Tag</a> for details. 
On these pages, there are lists of tuning parameters that can
potentially be optimized.  
<a href="custom_models.html">User-defined models</a> can also be created.
</p>

<p>The first step in tuning the model (line
1 in the algorithm above is to choose a set of
parameters to evaluate. For example, if fitting a Partial Least Squares
(PLS) model, the number of PLS components to evaluate must be specified. 
</p>

<p><br><img width = 629 height =234 src="TrainAlgo.png"><br><br></p>

<p>
Once the model and tuning parameter values have been defined, the type
of resampling should be also be specified. Currently, <i>k</i>-fold
cross-validation (once or repeated), 
leave-one-out cross-validation and bootstrap
(simple estimation or the 632 rule) 
resampling methods can be used by <span class="mx funCall">train</span>. After resampling,
the process produces a profile of performance measures is available to
guide the user as to which tuning parameter values should be
chosen. By default, the function automatically chooses the tuning
parameters associated with the best value, although different
algorithms can be used (see details below). 
</p>

<div id="example"></div> 
<h2>An Example</h2>

<p>The Sonar data are available in the <a href="http://cran.r-project.org/web/packages/mlbench/index.html"><strong>mlbench</strong></a> package. Here, we load the data:

<!--begin.rcode train_sonar1
library(mlbench)
data(Sonar)
str(Sonar[, 1:10])
    end.rcode-->

<p>
The function <span class="mx funCall">createDataPartition</span> can be used to create a stratified random sample of the data into training and test sets:
</p>
<!--begin.rcode train_sonar2
library(caret)
set.seed(998)
inTraining <- createDataPartition(Sonar$Class, p = .75, list = FALSE)
training <- Sonar[ inTraining,]
testing  <- Sonar[-inTraining,]
    end.rcode-->

<p>We will use these data illustrate functionality on this (and other) pages.
</p>

<div id="tune"></div> 
<h1>Basic Parameter Tuning</h1>

<p>
By default, simple bootstrap resampling is used for line
3 in the algorithm above. Others are availible, such
as repeated <i>K</i>-fold cross-validation, leave-one-out etc. The function
<span class="mx funCall">trainControl</span> can be used to specifiy the type of resampling:
</p>
<!--begin.rcode train_control,tidy=FALSE
fitControl <- trainControl(## 10-fold CV
                           method = "repeatedcv",
                           number = 10,
                           ## repeated ten times
                           repeats = 10)
    end.rcode-->
<p>
More information about <span class="mx funCall">trainControl</span> is given in  <a href="#custom">a section below</a>.
</p>

<p>The first two arguments to <span class="mx funCall">train</span> are the predictor and
outcome data objects, respectively. The third argument,
<span class="mx arg">method</span>, specifies the type of model (see <a href="modelList.html"><tt>train</tt> Model List</a> 
or <a href="bytag.html"><tt>train</tt> Models By Tag</a>). 
To illustrate, we will fit a boosted tree model via the <a href="http://cran.r-project.org/web/packages/gbm/index.html"><strong>gbm</strong></a>
package. The basic syntax for fitting this model using repeated
cross-validation is shown below:
</p>


<!--begin.rcode train_gbm1,cache=TRUE,tidy=FALSE
set.seed(825)
gbmFit1 <- train(Class ~ ., data = training, 
                 method = "gbm", 
                 trControl = fitControl,
                 ## This last option is actually one
                 ## for gbm() that passes through
                 verbose = FALSE)
gbmFit1
    end.rcode-->

<p>
For a gradient boosting machine (GBM) model, there are three main
tuning parameters:
  <ul>
  <li> number of iterations, i.e. trees,  (called <span class="mx arg">n.trees</span> in the
  <span class="mx funCall">gbm</span> function)</li>
  <li> complexity of the tree, called <span class="mx arg">interaction.depth</span></li>
  <li> learning rate: how quickly the algorithm adapts, called
  <span class="mx arg">shrinkage</span> </li>
  <li> the minimum number of training set samples in a node to commence splitting (<span class="mx arg">n.minobsinnode</span>)</li>
  </ul>
</p>  
<p>
The default values tested for this model are shown in the first two
columns (<span class="mx arg">shrinkage</span> and <span class="mx arg">n.minobsinnode</span> are not shown beause the grid set of
candidate models all use a single value for these tuning parameters). 
The column labeled "<code>Accuracy</code>" is the overall agreement rate
averaged over cross-validation iterations. The agreement standard
deviation is also calculated from the cross-validation results. The
column "<code>Kappa</code>" is Cohen's (unweighted) Kappa statistic
averaged across the resampling results. <span class="mx funCall">train</span> works with
specific models (see <a href="modelList.html"><tt>train</tt> Model List</a> 
or <a href="bytag.html"><tt>train</tt> Models By Tag</a>). 
For these models, <span class="mx funCall">train</span> can automatically
create a grid of tuning parameters. By default, if <i>p</i> is the number
of tuning parameters, the grid size is <i>3^p</i>. As another example, regularized
discriminant analysis (RDA) models have two  parameters
(<span class="mx arg">gamma</span> and <span class="mx arg">lambda</span>), both of which lie on [0,
1]. The default training grid would produce nine combinations in this
two-dimensional space. 
</p>

<p>There are several <a href="notes.html">notes</a> regarding specific model behaviors 
for <span class="mx funCall">train</span>. There is additional functionality in <span class="mx funCall">train</span> that is described in the next section.
</p>

<div id="repro"></div> 
<h1>Notes on Reproducibility</h1>

<p>
Many models utilize random numbers during the phase where parameters are estimated. Also, the resampling indices are chosen using random numbers. There are two main ways to control the randomness in order to assure reproducible results. 
</p>

<ul>
  <li>
  There are two approaches to ensuring that the same <em>resamples</em> are used between calls to <span class="mx funCall">train</span>. The first is to use <span class="mx funCall">set.seed</span> just prior to calling <span class="mx funCall">train</span>. The first use of random numbers is to create the resampling information. Alternatively, if you would like to use specific splits of the data, the <span class="mx arg">index</span> argument of the <span class="mx funCall">trainControl</span> function can be used. This is briefly discussed below.  
  </li>
  <li>
  When the models are created <em>inside of resampling</em>, the seeds can also be set. While setting the seed prior to calling <span class="mx funCall">train</span> may guarantee that the same random numbers are used, this is unlikely to be the case when <a href="parallel.html">parallel processing</a> is used (depending which technology is utilized). To set the model fitting seeds, <span class="mx funCall">trainControl</span> has an additional argument called <span class="mx arg">seeds</span> that can be used. The value for this argument is a list of integer vectors that are used as seeds. The help page for <span class="mx funCall">trainControl</span> describes the appropriate format for this option. 
  </li>
</ul>
<p>
How random numbers are used is highly dependent on the package author. There are rare cases where the underlying model function does not control the random number seed, especially if the computations are conducted in C code. Also, please note that <a href="https://github.com/topepo/caret/issues/452">some packages load random numbers when loaded (directly or via namespace)</a> and this may effect reproducibility.  
</p>

<div id="custom"></div> 
<h1>Customizing the Tuning Process</h1>

<p>
There are a few ways to customize the process of selecting
tuning/complexity parameters and building the final model.
</p>

<div id="preproc"></div> 
<h2>Pre-Processing Options</h2>

<p>
As previously mentioned,<span class="mx funCall">train</span> can pre-process the data in
various ways prior to model fitting. The function
<span class="mx funCall">preProcess</span> is automatically used. This function can be used
for centering and scaling, imputation (see details below),
applying the spatial sign transformation and feature extraction via
principal component analysis or independent component
analysis. 
</p>
<p>
To specify what pre-processing should occur, the <span class="mx funCall">train</span> function has an argument called <span class="mx arg">preProcess</span>. This argument takes a character string of methods that would normally be passed to the <span class="mx arg">method</span> argument of the <a href="preprocess.html"><span class="mx funCall">preProcess</span> function</a>. Additional options to the <span class="mx funCall">preProcess</span> function can be passed
via the <span class="mx funCall">trainControl</span> function. 
</p>
<p>
These processing steps would be applied during any predictions
generated using <span class="mx funCall">predict.train</span>, <span class="mx funCall">extractPrediction</span> or
<span class="mx funCall">extractProbs</span> (see details later in this document). The
pre-processing would <b>not</b> be applied to predictions that
directly use the <code>object$finalModel</code> object.  
</p>
<p>
For imputation, there are three methods currently implemented:
</p>

  <ul>
  <li> <i>k</i>-nearest neighbors takes a sample with missing values and
  finds the <i>k</i> closest samples in the training set. The average of
  the <i>k</i> training set values for that predictor are used as a
  substitute for the original data. When calculating the distances to
  the training set samples, the predictors used in the calculation are
  the ones with no missing values for that sample and no missing
  values in the training set.</li>
  <li> another approach is to fit a bagged tree model for each
  predictor using the training set samples. This is usually a fairly
  accurate model and can handle missing values. When a predictor for a
  sample requires imputation, the values for the other predictors are
  fed through the bagged tree and the prediction is used as the new
  value. This model can have significant computational cost.</li>
  <li> the median of the predictor's training set values can be used
  to estimate the missing data.</li>  
  </ul>
<p>
If there are missing values in the training set, PCA and ICA models
only use complete samples.
</p>

<div id="grids"></div> 
<h2>Alternate Tuning Grids</h2>
<p>
The tuning parameter grid can be specified by the user. The argument
<span class="mx arg">tuneGrid</span> can take a data frame with columns for each tuning parameter. The column names should be the same as the fitting function's arguments. For the previously mentioned RDA example, the names would be
<code>gamma</code> and <code>lambda</code>. <span class="mx funCall">train</span> will tune the
model over each combination of values in the rows. 
</p>
<p>
For the boosted tree model, we can fix the learning rate and evaluate more than three values of
<span class="mx arg">n.trees</span>: 
</p>

<!--begin.rcode train_gbm2,cache=TRUE,tidy=FALSE,results='hide'
gbmGrid <-  expand.grid(interaction.depth = c(1, 5, 9), 
                        n.trees = (1:30)*50, 
                        shrinkage = 0.1,
                        n.minobsinnode = 20)
                        
nrow(gbmGrid)

set.seed(825)
gbmFit2 <- train(Class ~ ., data = training, 
                 method = "gbm", 
                 trControl = fitControl, 
                 verbose = FALSE, 
                 ## Now specify the exact models 
                 ## to evaluate:
                 tuneGrid = gbmGrid)
gbmFit2
    end.rcode-->
<!--begin.rcode train_gbm2_print,echo=FALSE
gbmOut <- capture.output(print(gbmFit2, digits = 2))
text2 <- c(gbmOut[1:18],
           "  :                   :        :         : ",
           gbmOut[(length(gbmOut)-10):length(gbmOut)])
cat(paste(text2, collapse = "\n"))
    end.rcode-->

<p>
If there are missing values in the training set, PCA and ICA models
only use complete samples.
</p>

<p>
Another option is to use a random sample of possible tuning parameter combinations, i.e. "random search"<a href="www.jmlr.org/papers/volume13/bergstra12a/bergstra12a.pdf">(pdf)</a>. This functionality is described on <a href="random.html">this page</a>. 
</p>

<p>
To use a random search, use the option <span class="mx arg">search</span><tt> = <!--rinline '"random"' --></tt> in the call to <span class="mx funCall">trainControl</span>. In this situation, the <span class="mx arg">tuneLength</span> parameter defines the total number of parameter combinations that will be evaluated. 
</p>


<div id="plots"></div> 
<h2>Plotting the Resampling Profile</h2>

<p>The <span class="mx funCall">plot</span> function can be used to examine the relationship between the 
estimates of performance and the tuning parameters. For example, a simple invokation of 
the function shows the results for the first performance measure:
</p>


<!--begin.rcode train_plot1,fig.width=7,fig.height=4  
trellis.par.set(caretTheme())
plot(gbmFit2)  
    end.rcode-->
    
<p>Other performance metrics can be shown using the <span class="mx arg">metric</span> option:


<!--begin.rcode train_plot2,fig.width=7,fig.height=4
trellis.par.set(caretTheme())
plot(gbmFit2, metric = "Kappa")
    end.rcode-->

<p>Other types of plot are also available. See <code>?plot.train</code> for more details.
The code below shows a heatmap of the results:
</p>


<!--begin.rcode train_plot3,tidy=FALSE,fig.width=7,fig.height=4
trellis.par.set(caretTheme())
plot(gbmFit2, metric = "Kappa", plotType = "level",
     scales = list(x = list(rot = 90)))
    end.rcode-->

    
<p>A <span class="mx funCall">ggplot</span> method can also be used:</p>

<!--begin.rcode train_ggplot1,fig.width=8.5,fig.height=4  
ggplot(gbmFit2)  
    end.rcode-->


<p>There are also plot functions that show more detailed representations of the 
resampled estimates. See <code>?xyplot.train</code> for more details.
</p>
<p>
From these plots, a different set of tuning parameters may be desired. To change the 
final values without starting the whole process again, the <span class="mx funCall">update.train</span>
can be used to refit the final model. See <code>?update.train</code>
</p>

<div id="control"></div> 
<h2>The <span class="mx funCall">trainControl</span> Function</h2>

<p>
The function <span class="mx funCall">trainControl</span> generates parameters that further
control how models are created, with possible values: 
</p>
<ul>
 <li> <span class="mx arg">method</span>: The resampling method: <tt><!--rinline '"boot"' --></tt>,  <tt><!--rinline '"cv"' --></tt>, <tt><!--rinline '"LOOCV"' --></tt>, <tt><!--rinline '"LGOCV"' --></tt>, <tt><!--rinline '"repeatedcv"' --></tt>, <tt><!--rinline '"timeslice"' --></tt>, <tt><!--rinline '"none"' --></tt> and <tt><!--rinline '"oob"' --></tt>. 
 The
  last value, out-of-bag estimates, can only be used by random
  forest, bagged trees, bagged earth, bagged flexible discriminant
  analysis, or conditional tree forest models. GBM models are not
  included (the <a href="http://cran.r-project.org/web/packages/gbm/index.html"><strong>gbm</strong></a> package maintainer has indicated that
  it would not be a good idea to choose tuning parameter values
  based on the model OOB error estimates with boosted trees). Also,
  for leave-one-out cross-validation, no uncertainty estimates
  are given for the resampled performance measures.  </li>
 <li> <span class="mx arg">number</span> and <span class="mx arg">repeats</span>: <span class="mx arg">number</span> controls with the 
  number of folds in <i>K</i>-fold cross-validation or number of
  resampling iterations for bootstrapping and leave-group-out
  cross-validation. <span class="mx arg">repeats</span> applied only to repeated  
  <i>K</i>-fold cross-validation. Suppose that 
  <span class="mx arg">method</span><tt> = <!--rinline '"repeatedcv"' --></tt>,
  <span class="mx arg">number</span><tt> = <!--rinline '10' --></tt> and
  <span class="mx arg">repeats</span><tt> = <!--rinline '3' --></tt>,then three separate
  10-fold cross-validations are used as the resampling scheme. </li>
 <li> <span class="mx arg">verboseIter</span>: A logical for printing a training log. </li>
 <li> <span class="mx arg">returnData</span>: A logical for saving the data into a slot
  called <code>trainingData</code>.  </li>
 <li> <span class="mx arg">p</span>: For leave-group out cross-validation: the training
  percentage  </li>
  <li>For <span class="mx arg">method</span><tt> = <!--rinline '"timeslice"' --></tt>, <span class="mx funCall">trainControl</span> has options <span class="mx arg">initialWindow</span>, <span class="mx arg">horizon</span> and <span class="mx arg">fixedWindow</span> that govern how <a href="splitting.html">cross-validation can be used for time series data. </a> 
  <li> <span class="mx arg">classProbs</span>: a logical value determining whether
  class probabilities should be computed for held-out samples
  during resample.  </li>
 <li> <span class="mx arg">index</span> and <span class="mx arg">indexOut</span>: optional lists with elements for each resampling
  iteration. Each list element is the sample rows used for training
  at that iteration or should be held-out. When these values are not specified,
  <span class="mx funCall">train</span> will generate them.  </li>
 <li> <span class="mx arg">summaryFunction</span>: a function to compute alternate
  performance summaries. </li> 
 <li> <span class="mx arg">selectionFunction</span>: a function to choose the optimal
  tuning parameters. 
  and examples.  </li>
 <li> <span class="mx arg">PCAthresh</span>, <span class="mx arg">ICAcomp</span> and <span class="mx arg">k</span>: these are
  all options to pass to the <span class="mx funCall">preProcess</span> function (when used). </li>
 <li> <span class="mx arg">returnResamp</span>: a character string containing one of
  the following values: <tt><!--rinline '"all"' --></tt>, <tt><!--rinline '"final"' --></tt> or
  <tt><!--rinline '"none"' --></tt>. This specifies how much of the resampled
  performance measures to save. </li>  
  <li> <span class="mx arg">allowParallel</span>: a logical that governs whether <span class="mx funCall">train</span> should <a href="parallel.html">use parallel processing (if availible). </a>  </li>
</ul>
<p>There are several other options not discussed here.<p/>

<div id="metrics"></div> 
<h2>Alternate Performance Metrics</h2>

<p>
The user can change the metric used to determine the best settings. By
default, RMSE and <i>R</i><sup>2</sup> are computed for regression while accuracy and
Kappa are computed for classification. Also by default, the parameter
values are chosen using RMSE and accuracy, respectively  for
regression and classification. The <span class="mx arg">metric</span> argument of the
<span class="mx funCall">train</span> function allows the user to control which the
optimality criterion is used. For example, in problems where there are
a low percentage of samples in one class, using <span class="mx arg">metric</span><tt> = <!--rinline '"Kappa"' --></tt> can improve quality of the final model. 
</p>
<p>
If none of these parameters are satisfactory, the user can also
compute custom performance metrics. The <span class="mx funCall">trainControl</span> function
has a argument called <span class="mx arg">summaryFunction</span> that specifies a
function for computing performance. The function should have these
arguments: 
</p>
<p>
<ul>
 <li> <span class="mx arg">data</span> is a reference for a data frame or matrix with
  columns called <code>obs</code> and <code>pred</code> for the observed and
  predicted outcome values (either numeric data for regression or
  character values for classification). Currently, class probabilities
  are not passed to the function. The values in data are the held-out
  predictions (and their associated reference values) for a single
  combination of tuning parameters. If the <span class="mx arg">classProbs</span>
  argument of the <span class="mx funCall">trainControl</span> object is set to
  <tt><!--rinline 'TRUE' --></tt>, additional columns in <code>data</code> will be present
  that contains the class probabilities. The names of these columns
  are the same as the class levels. Also, if  <span class="mx arg">weights</span> were specified in the call to <span class="mx funCall">train</span>, a column called <code>weights</code> will also be in the data set. </li>
 <li> <span class="mx arg">lev</span> is a character string that has the outcome factor
  levels taken from the training data. For regression, a value of
  <tt><!--rinline 'NULL' --></tt> is passed into the function. </li>
 <li> <span class="mx arg">model</span> is a character string for the model being used
  (i.e. the value passed to the <span class="mx arg">method</span> argument of
  <span class="mx funCall">train</span>). </li>
</ul>
</p>
<p>
The output to the function should be a vector of numeric summary
metrics with non-null names.  By default, <span class="mx funCall">train</span> evaluate classification models in terms of
the predicted classes. Optionally, class probabilities can also be
used to measure performance. To obtain predicted class probabilities
within the resampling process, the argument <span class="mx arg">classProbs</span> in
<span class="mx funCall">trainControl</span> must be set to <tt><!--rinline 'TRUE' --></tt>. This merges
columns of probabilities into the predictions generated from each
resample (there is a column per class and the column names are the
class names).
</p>
<p>
As shown in the last section, custom functions can be used to
calculate performance scores that are averaged over the
resamples. Another built-in function, <span class="mx funCall">twoClassSummary</span>, will
compute the sensitivity, specificity and area under the ROC curve:
</p>
<!--begin.rcode train_summary1
head(twoClassSummary)
    end.rcode-->

<p>
To rebuild the boosted tree model using this criterion, we
can see the relationship between the tuning parameters and
the area under the ROC curve using the following code: 
</p>
               
<!--begin.rcode train_summary2,cache=TRUE,tidy=FALSE,results='hide'
fitControl <- trainControl(method = "repeatedcv",
                           number = 10,
                           repeats = 10,
                           ## Estimate class probabilities
                           classProbs = TRUE,
                           ## Evaluate performance using 
                           ## the following function
                           summaryFunction = twoClassSummary)

set.seed(825)
gbmFit3 <- train(Class ~ ., data = training, 
                 method = "gbm", 
                 trControl = fitControl, 
                 verbose = FALSE, 
                 tuneGrid = gbmGrid,
                 ## Specify which metric to optimize
                 metric = "ROC")
gbmFit3
    end.rcode-->
<!--begin.rcode train_gbm3_print,echo=FALSE
gbmOut <- capture.output(print(gbmFit3, digits = 2))
text3 <- c(gbmOut[1:18],
           "  :                   :        :     :      :    ",
           gbmOut[(length(gbmOut)-10):length(gbmOut)])
cat(paste(text3, collapse = "\n"))
    end.rcode-->    
<!--begin.rcode train_pred3,cache=TRUE,echo=FALSE,results='hide'
testPred <- predict(gbmFit3, testing) ## to be used in 'other'html'
    end.rcode--> 
<p>
In this case, the average area under the ROC curve associated with the
optimal tuning parameters was
<!--rinline I(round(caret:::getTrainPerf(gbmFit3)[1,"TrainROC"], 3)) --> across
the <!--rinline I(length(gbmFit3$control$index)) --> resamples.
</p>

<div id="final"></div> 
<h2>Choosing the Final Model</h2>

<p>
Another method for customizing the tuning process is to modify the
algorithm that is used to select the "best" parameter values, given
the performance numbers. By default, the <span class="mx funCall">train</span> function
chooses the model with the largest performance value (or smallest, for
mean squared error in regression models). Other schemes for selecting
model can be used.  <a href="http://books.google.com/books/about/Classification_and_Regression_Trees.html?id=JwQx-WOmSyQC">Breiman et al (1984)</a> suggested the "one standard
error rule" for simple tree-based models. In this case, the model
with the best performance value is identified and, using resampling,
we can estimate the standard error of performance. The final model
used was the simplest model within one standard error of the
(empirically) best model. With simple trees this makes sense, since
these models will start to over-fit as they become more and more
specific to the training data. 
</p>
<p>
<span class="mx funCall">train</span> allows the user to specify alternate rules for
selecting the final model. The argument <span class="mx arg">selectionFunction</span>
can be used to supply a function to algorithmically determine the
final model. There are three existing functions in the package:
<span class="mx arg">best</span> is chooses the largest/smallest value, <span class="mx funCall">oneSE</span>
attempts to capture the spirit of <a href="http://books.google.com/books/about/Classification_and_Regression_Trees.html?id=JwQx-WOmSyQC">Breiman et al (1984)</a> and
<span class="mx funCall">tolerance</span> selects the least complex model within some percent
tolerance of the best value. See <code>?best</code> for more details. 
</p>
<p>
User-defined functions can be used, as long as they have the
following arguments: 
<ul>
  <li><span class="mx arg">x</span> is a data frame containing the tune parameters and
    their associated performance metrics. Each row corresponds to a
    different tuning parameter combination. </li> 
  <li><span class="mx arg">metric</span> a character string indicating which
      performance metric should be optimized (this is passed in
      directly from the <span class="mx arg">metric</span> argument of <span class="mx funCall">train</span>.</li>  
  <li><span class="mx arg">maximize</span> is a single logical value indicating
        whether larger values of the performance metric are better
        (this is also directly passed from the call to
        <span class="mx funCall">train</span>). </li>
</ul>
<p>
  The function should output a single integer indicating which row in
  <code>x</code> is chosen. 
</p>
<!--begin.rcode bestGBM,echo=FALSE
printSelected <- function(x) {
    tmp <- x$bestTune
    tmp <- paste(names(tmp), "=", tmp)
    paste(tmp, collapse = ", ")
  }
    end.rcode-->   
  
<p>  
As an example, if we chose the previous boosted tree model on the basis of
overall accuracy, we would choose:
<!--rinline I(printSelected(gbmFit3)) -->. However, the scale in this plots is
fairly tight, with accuracy values ranging from
<!--rinline I(round(min(gbmFit3$results$ROC), 3)) --> to
<!--rinline I(round(max(gbmFit3$results$ROC), 3)) -->. A less complex model
(e.g. fewer, more shallow trees) might also yield acceptable
accuracy. 
</p>
<p>
The tolerance function could be used to find a less complex
model based on (<i>x</i>-<i>x</i><sub>best</sub>)/<i>x</i><sub>best</sub>x 100, which is the percent
difference. For example, to select parameter values based on a 2% loss of performance: 
</p>
<!--begin.rcode train_tolerance
whichTwoPct <- tolerance(gbmFit3$results, metric = "ROC", 
                         tol = 2, maximize = TRUE)  
cat("best model within 2 pct of best:\n")
gbmFit3$results[whichTwoPct,1:6]
    end.rcode-->
    
<p>
This indicates that we can get a less complex model with an area under the ROC curve 
of <!--rinline I(round(gbmFit3$results[whichTwoPct,"ROC"], 3)) --> (compared
to the "pick the best" value of 
 <!--rinline I(round(getTrainPerf(gbmFit3)[1,"TrainROC"], 3)) -->). 
</p>
<p>
The main issue with these functions is related to ordering the models
from simplest to complex. In some cases, this is easy (e.g. simple
trees, partial least squares), but in cases such as this model, the ordering of
models is subjective. For example, is a boosted tree model using 100
iterations and a tree depth of 2 more complex than one with 50
iterations and a depth of 8? The package makes some choices regarding
the orderings. In the case of boosted trees, the package assumes that
increasing the number of iterations adds complexity at a faster rate
than increasing the tree depth, so models are ordered on the number of
iterations then ordered with depth. See <code>?best</code> for more
examples for specific models. 
</p>

<div id="pred"></div> 
<h1>Extracting Predictions and Class Probabilities</h1>

<p>
As previously mentioned, objects produced by the <span class="mx funCall">train</span>
function contain the "optimized" model in the <code>finalModel</code>
sub-object. Predictions can be made from these objects as usual. In
some cases, such as <span class="mx funCall">pls</span> or <span class="mx funCall">gbm</span> objects, additional
parameters from the optimized fit may need to be specified. In these
cases, the <span class="mx funCall">train</span> objects uses the results of the parameter
optimization to predict new samples. For example, if predictions were create
using <span class="mx funCall">predict.gbm</span>, the user would have to specify the number of trees 
directly (there is no default). Also, for binary classification, the predictions
from this function take the form of the probability of one of the classes, so 
extra steps are required to convert this to a factor vector. <span class="mx funCall">predict.train</span>
automatically handles these details for this (and for other models).
</p>
<p>
Also, there are very few standard syntaxes for model predictions in R. For example, 
to get class probabilities, many <span class="mx funCall">predict</span> methods have an argument called 
<span class="mx arg">type</span> that is used to specify whether the classes or probabilities should 
be generated. Different packages use different values of <span class="mx arg">type</span>, such as 
<tt><!--rinline '"prob"' --></tt>, <tt><!--rinline '"posterior"' --></tt>, <tt><!--rinline '"response"' --></tt>, <tt><!--rinline '"probability"' --></tt> or <tt><!--rinline '"raw"' --></tt>. In other cases, completely 
different syntax is used.
</p>
<p>For <span class="mx funCall">predict.train</span>, the type options are standardized to be <tt><!--rinline '"class"' --></tt> and <tt><!--rinline '"prob"' --></tt> (the underlying code matches these to the appropriate choices for each model. For example:
</p>
<!--begin.rcode train_pred1
predict(gbmFit3, newdata = head(testing))

predict(gbmFit3, newdata = head(testing), type = "prob")
    end.rcode-->

<div id="resamp"></div> 
<h1>Exploring and Comparing Resampling Distributions</h1>

<div id="within"></div> 
<h2>Within-Model</h2>
<p>
There are several <a href="http://cran.r-project.org/web/packages/lattice/index.html"><strong>lattice</strong></a> functions than can be used to explore
relationships between tuning parameters and the resampling results for
a specific model:
<ul>
 <li><span class="mx funCall">xyplot</span> and <span class="mx funCall">stripplot</span> can be used to plot resampling statistics
  against (numeric) tuning parameters.</li>
 <li><span class="mx funCall">histogram</span> and <span class="mx funCall">densityplot</span> can also be used
      to look at distributions of the tuning parameters across tuning parameters.</li>
</ul>
</p>
<p>
For example, the following statements create a density plot: 
</p>


<!--begin.rcode 4,echo=FALSE,fig.width=7,fig.height=4
trellis.par.set(caretTheme())
densityplot(gbmFit3, pch = "|")
    end.rcode-->

<p>Note that if you are interested in plotting the resampling results across multiple tuning parameters, the option <tt><!--rinline 'resamples = "all"' --></tt> should be used in the control object.</p>

<div id="between"></div> 
<h2>Between-Models</h2>

<p>
The <a href="http://cran.r-project.org/web/packages/caret/index.html"><strong>caret</strong></a> package also includes functions to characterize the differences
between models (generated using <span class="mx funCall">train</span>, <span class="mx funCall">sbf</span> or
<span class="mx funCall">rfe</span>) via their resampling distributions. These functions are
based on the work of <a href="http://www.stat.uni-muenchen.de/~leisch/papers/Hothorn+Leisch+Zeileis-2005.pdf">Hothorn et al. (2005)</a> and <a href="http://epub.ub.uni-muenchen.de/10604/1/tr56.pdf">Eugster et al (2008)</a>.
</p>
<p>
First, a support vector machine model is fit to the Sonar data. The data are centered and scaled using the <span class="mx arg">preProc</span> argument. Note that the same random number seed is set prior to the model that is identical to the seed used for the boosted tree model. This ensures that the same resampling sets are used, which will come in handy when we compare the resampling profiles between models. </p>
</p>
<!--begin.rcode train_svmFit,cache=TRUE,tidy=FALSE
set.seed(825)
svmFit <- train(Class ~ ., data = training, 
                 method = "svmRadial", 
                 trControl = fitControl, 
                 preProc = c("center", "scale"),
                 tuneLength = 8,
                 metric = "ROC")
svmFit                 
    end.rcode-->
<p>
Also, a regularized discriminant analysis model was fit.  
</p>
<!--begin.rcode train_rdaFit,cache=TRUE,tidy=FALSE
set.seed(825)
rdaFit <- train(Class ~ ., data = training, 
                 method = "rda", 
                 trControl = fitControl, 
                 tuneLength = 4,
                 metric = "ROC")
rdaFit                 
    end.rcode-->
<p>
Given these models, can we make statistical statements about their
performance differences? To do this, we first collect the resampling
results using <span class="mx funCall">resamples</span>. 
</p>
<!--begin.rcode train_resamps1,tidy=FALSE
resamps <- resamples(list(GBM = gbmFit3,
                          SVM = svmFit,
                          RDA = rdaFit))
resamps

summary(resamps)
    end.rcode-->

<p>Note that, in this case, the option <tt><!--rinline 'resamples = "final"' --></tt> should be used in the control objects.</p>

<p>
There are several <pkg>lattice</pkg> plot methods that can be used to visualize
the resampling distributions: density plots, box-whisker plots,
scatterplot matrices and scatterplots of summary statistics. For example:
</p>

<!--begin.rcode train_resample_box,fig.width=9,fig.height=4
trellis.par.set(theme1)
bwplot(resamps, layout = c(3, 1))
    end.rcode-->

<!--begin.rcode train_resample_ci,fig.height=6,fig.width=6
trellis.par.set(caretTheme())
dotplot(resamps, metric = "ROC")
    end.rcode-->

<!--begin.rcode train_resample_ba,fig.height=6,fig.width=6
trellis.par.set(theme1)
xyplot(resamps, what = "BlandAltman")
    end.rcode-->


<!--begin.rcode train_resample_scatmat
splom(resamps)
    end.rcode-->

<p> Other visualizations are availible in <span class="mx funCall">densityplot.resamples</span> and <span class="mx funCall">parallel.resamples</span></p>


<p>
Since models are fit on the same versions of the training data, it
makes sense to make inferences on the differences between models. In
this way we reduce the within-resample correlation that may exist. We
can compute the differences, then use a simple <i>t</i>-test to evaluate
the null hypothesis that there is no difference between models. 
</p>
<!--begin.rcode train_resamps3
difValues <- diff(resamps)
difValues

summary(difValues)
    end.rcode-->


<!--begin.rcode train_diff_box,fig.width=9,fig.height=4
trellis.par.set(theme1)
bwplot(difValues, layout = c(3, 1))
    end.rcode-->


<!--begin.rcode train_diff_ci,fig.height=6,fig.width=6
trellis.par.set(caretTheme())
dotplot(difValues)
    end.rcode-->

<div id="notune"></div> 
<h1>Fitting Models Without Parameter Tuning</h1>

<p>In cases where the model tuning values are known, <span class="mx funCall">train</span> can be used to fit the model to the entire training set without any resampling or parameter tuning. Using the <tt><!--rinline 'method = "none"' --></tt> option in <span class="mx funCall">trainControl</span> can be used. For example:
</p>
<!--begin.rcode train_noresamp,cache=TRUE,tidy=FALSE
fitControl <- trainControl(method = "none", classProbs = TRUE)

set.seed(825)
gbmFit4 <- train(Class ~ ., data = training, 
                 method = "gbm", 
                 trControl = fitControl, 
                 verbose = FALSE, 
                 ## Only a single model can be passed to the
                 ## function when no resampling is used:
                 tuneGrid = data.frame(interaction.depth = 4,
                                       n.trees = 100,
                                       shrinkage = .1,
                                       n.minobsinnode = 20),
                 metric = "ROC")
gbmFit4
    end.rcode-->

<p>
Note that <span class="mx funCall">plot.train</span>, <span class="mx funCall">resamples</span>, <span class="mx funCall">confusionMatrix.train</span> and several other functions will not work with this object but <span class="mx funCall">predict.train</span> and others will:
</p>
<!--begin.rcode train_noresamppred1
predict(gbmFit4, newdata = head(testing))
predict(gbmFit4, newdata = head(testing), type = "prob")
    end.rcode-->

<div style="clear: both;">&nbsp;</div>
  </div>
  <!-- end #content -->
<div id="sidebar">
<ul>
  <li>
    <h2>General Topics</h2>
    <ul>
      <li><a href="index.html">Front Page</a></li>
      <li><a href="visualizations.html">Visualizations</a></li>
      <li><a href="preprocess.html">Pre-Processing</a><li>
      <li><a href="splitting.html">Data Splitting</a></li>
      <li><a href="varimp.html">Variable Importance</a></li>
      <li><a href="other.html">Model Performance</a></li>
      <li><a href="parallel.html">Parallel Processing</a></li>
    </ul>
    <h2>Model Training and Tuning</h2>
    <ul>
      <li><a href="training.html">Basic Syntax</a></li>
      <li><a href="modelList.html">Sortable Model List</a></li>
      <li><a href="bytag.html">Models By Tag</a></li>
      <li><a href="similarity.html">Models By Similarity</a></li>
      <li><a href="custom_models.html">Using Custom Models</a></li>
      <li><a href="sampling.html">Sampling for Class Imbalances</a></li> 
      <li><a href="random.html">Random Search</a></li> 
      <li><a href="adaptive.html">Adaptive Resampling</a></li> 
    </ul>
    <h2>Feature Selection</h2>
    <ul>
      <li><a href="featureselection.html">Overview</a>
      <li><a href="rfe.html">RFE</a></li>
      <li><a href="filters.html">Filters</a></li>
      <li><a href="GA.html">GA</a></li>
      <li><a href="SA.html">SA</a></li>
    </ul>  
  </li>
</ul>
</div>
<!-- end #sidebar -->
<div style="clear: both;">&nbsp;</div>
  </div>
  <div class="container"><img src="images/img03.png" width="1000" height="40" alt="" /></div>
  <!-- end #page -->
</div>
<!--begin.rcode echo = FALSE
knit_hooks$set(inline = hook_inline)    
    end.rcode--> 
  <div id="footer-content"></div>
  <div id="footer">
  <p>Created on <!--rinline I(session) -->.</p>
  </div>
  <!-- end #footer -->
</body>
  </html>
